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ABSTRACT 


A fast  search  algorithm  for  the  solution  of  nonlinear  mathemati- 
cal programming  optimization  problems  is  presented  in  this  Report . A 
gradient  search  procedure  is  combined  with  a "Boundary  Tracking "(BT) 
method  using  the  feasible  direction  finding  method  of  Zoutendijk  for 
generating  a feasible  starting  direction  along  the  feasible-infeasible 
boundary. 

The  algorithm  is  applied  to  the  minimum  weight  design  of  submersi- 
ble, circular,  cylindrical  shells  reinforced  by  equally  spaced  "T"  type 
frames.  This  problem  had  produced  algorithm  failure  in  two  earlier 
studies  and  was  only  recently  solved  by  the  Direct  Search-Feasible 
Direction  Algorithm  (DSFD)  which  was  shown  by  recent  comparison  studies 
to  be  among  the  fastest  and  most  reliable  mathematical  programming 
methods  available.  The  BT  procedure  was  found  to  be  substantially 
faster  than  DSFD,  producing  a solution  with  about  one-eigth  the  effort 
required  by  DSFD. 

In  a general  comparison  study  a code  based  on  the  BT  algorithm 
was  compared  with  twenty  other  codes  representing  most  of  the  popular 
numerical  optimization  methods  on  ten  test  problems.  These  problems 
are  such  that  majority  of  the  codes  tested  failed  to  solve  more  than 
half  of  them.  The  new  code  proved  superior  to  all  others  in  overall 
generality  and  efficiency.  It  solved  all  problems  and  was  the  fastest 
code  on  the  constrained  problems. 
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NOMENCLATURE 
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A = cross  sectional  area  of  the  frame 

A.  = number  of  active  behavior  constraints 

w • 

Bk(x)  = behavior  function 

b = web  thickness 

b'  = flange  thickness 

c = distance  from  midplane  of  the  shell  to  the  surface  of 

the  frame 

C-]  = an  arbitrary  large  constant 

C2  = maximum  flange  to  plating  thickness  ratio 

d = distance  from  the  midplane  of  the  hull  plating  to  the 

neutral  axis  of  the  frame 

E = tensile  modulus 

EIf  = bending  stiffness  of  a frame 

e = eccentricity 

e-j,  e2>  e3  = arbitrary  small  positive  constraints 
f,  f(x)  = objective  function 

f . T = efficiency  criteria 

a a 

G = shear  modulus 

GJ^  = torsional  stiffness  of  a frame 

g(x),  g.(x)  = constraint  function  and  the  jth  constraint  function 

J 

g*(x)  = the  smallest  proper  constraint 

g(t)  „ = the  constraint  g*(x)  as  a function  of  a single  variable 

t 

Hq,  Hm  = frame  deflection  parameters  [see  equations  3.11  - 3.12] 
h = web  height 


I 


= number  of  variables 


J 

K 

kj 

L' 

Lk 

Ls 

m 

Na 

Nr  N2 
n 

na 

P 

P, 


PCg(n’m) 

Pcs(n,m) 

p*  , P* 
Hcg’  Kcs 

% 

R 

R R 
min’  max 

Ri’  Ri 

Sj 

t 

t 


ap 


number  of  constraint  equations 

number  of  proper  constraints 

arbitrary  small  positive  number 

distance  between  frames 

lower  limit  on  behavior  function 

length  of  the  shell 

axial  wave  number 

numerical  success  rating 

arbitrary  selected  numbers 

circumferential  wave  number 

number  of  problems  solved  by  code  "a" 

hydrostatic  pressure 

set  of  proper  constraints 

progress  towards  solution  (section  4.2) 

gross  collapse  pressure 

collapse  pressure  of  a shell  panel 

minimum  collapse  pressures 

total  radial  load 

radius  of  the  shell 

minimum  and  maximum  radius  of  the  shell 

active  lower  arid  upper  regional  constraints 

factor  of  safety 

skin  thickness 

normalized  cpu  time 

best  movement  direction  vector 


displacement  volume  of  the  cylinder 

volume  of  the  frame 

volume  of  plating 

volume  of  frame  web 

weight  of  the  hull 

deflection  parameter 

design  variable  vector  and  its  components  respectively 

step  size 

minimum  step  size 

frame  deflection  parameter 

specific  weight  of  material 

specific  weight  of  immersion  fluid 

constraint  activity  limit 

poison's  ratio 

frame  bending  stress 

compressive  hoop  stress 

allowable  frame  stress  . 

allowable  plating  stress 

axial  stress 

maximum  frame  stress 

hoop  stress 

magnitude  of  arbitrary  vector  function  $ 
the  gradient  of  arbitrary  vector  function  4) 
transpose  of  arbitrary  vector  function  <fi 
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SUPERSCRIPTS 

b = base  number 

l = comparison  value 

L = Lower  limit 

'Hi 

n = value  at  the  n direction  finding  problem 

r = number  of  redesign  cycles 

s = number  of  iterations  for  locating  the  IF  boundary 


U 


= upper  limit 
initial 
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INTRODUCTION 

1.1  Optimization  in  Design 

Optimal  design  problems  may  be  treated  by  many  different  methods, 
for  example,  ordinary  and  variational  calculus,  mathematical  program- 
ming  (MP),  and  some  special  techniques,  such  as  the  fully  stressed 
concept  used  in  structural  design  [1-10]\ 

Of  the  above  mentioned  methods,  MP  procedures  [11]  seem  to  be 
the  most  flexible  methods  available  for  optimal  design  synthesis 
since  they  treat  the  broadest  range  of  engineering  problems  and  are 
easily  adaptable. 

The  concept  of  design  optimization  requires  that  a quantity  to 
be  minimized  or  maximized  be  designated  and  that  this  quantity  be 
expressed  as  a function  of  the  design  variables.  This  function  is 
called  the  "merit"  or  "objective"  function,  and  can  be  written  in 
the  form 

f = f (x.)  1 = 1,2, ...I  (1-1) 

where  x-  are  the  design  variables  and  I is  the  number  of  design  vari- 
ables. Generally,  the  design  variables  are  not  free  to  take  on  any 
value,  but  are  subject  to  constraints  governing  their  range  or  the 
behavior  of  the  design. 
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Numbers  in  brackets  designate  references  at  the  end  of  the  thesis, 


, 
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When  limits  on  behavior  are  specified,  a quantitative  measure 
of  the  behavior  as  a function  of  the  variables,  called  the  "behavior 
function",  B(xl<)  is  required.  Thus,  if  K behavior  functions  are 
specified,  one  can  write  K equations  of  the  form 
Lk  Uk  k = 1,  2, . . .m 

Bk  < Uk  k = m + 1,  m + 2,...n  (1.2) 

l-k  < Bk  k = n + 1,  n + 2,...p 

where  Lk  is  the  lower  bound  on  Bk,  and  Uk  the  upper  limit.  !_k  and 

Uk  may  be  functions  of  x...  Behavior  constraints  representing  the 
above  behavior  function  limits  can  be  written  as  follows 

9j  = Uj  - Bj  > ° j = 1,  2,...m 

9J  ’ Bj-m  - Lj-m  5 0 j * m + 1 , m + 2,.. .2m  (1.3) 

gj  = Uj-2m  ' Bj-2m  ^ 0 j = 2m  + 1 , 2m  + 2, . . .2m  + n 

9j  = Bj-2m+n  - Lj-2m+n  * 0 j = 2m  + n + l,...2m  ...p 

All  the  above  equations  may  be  written  as 

9j  * 0 j = 1,  2....J  (1.4) 

J = 2m  + n + p 

Regional  limits  are  often  imposed  by  manufacturing  or  other 
considerations.  These  are  of  the  form 
L U 

x.  4 x.  ^ xi  (1 .5) 

x!r  and  x'j1  are  constants  and  represent  the  min-imum  and  maximum  values 
of  x..  respectively.  Here  regional  constraints  are  distinguished  from 
behavior  constraints  since  any  constraint  of  the  form 
b 4 xi  < B 


(1.6) 
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where  A and  B are  constants,  may  be  treated  more  simply  than  the 
general  form.  Not  all  the  x..  need  be  subject  to  such  limits.  A de- 
tailed discussion  of  regional  constraints  can  be  found  in  Reference 
[9]. 

Some  variables  may  also  be  restricted  to  certain  discrete  values. 

In  structural  design  these  values  can  represent  material  properties, 
geometric  available  shapes  and  sizes,  or  thickness  gages. 

The  concept  of  a merit  surface  is  important  to  the  understanding 
of  the  redesign  process.  The  function  fCx^ ) can  be  considered  as  a sur- 
face in  I + 1 dimensional  space  where  the  coordinate  axes  are  the  x^  and 
y where  the  value  of  y is  given  by  f(x^).  The  constraints  delineate  the 
region  of  interest  within  which  the  optimum  is  to  be  found.  Points 
which  satisfy  the  constraint  equations  are  called  "acceptable"  or 
"feasible"  points.  All  other  points  are  called  "unacceptable"  or 
"infeasible".  The  set  of  all  feasible  points  constitute  the  "feasible 
region"  and  all  infeasible  points  define  the  "infeasible  region".  The 
surface  between  these  regions  is  called  the  Infeasible-Feasible  (IF) 
boundary.  Those  portions  of  the  IF  boundary  where  at  least  one  be- 
havior constraint  is  zero  are  called  the  "behavior"  boundary.  Points 
away  from  the  constraint  surface  are  called  "free",  and  those  on 
these  surfaces  are  called  "bound"  points.  The  merit  surface  is  ex- 
plored to  find  the  highest  point  in  the  acceptable  region.  The  algo- 
rithms commonly  employed  in  such  problems  usually  start  the  exploration 
from  a free  acceptable  point.  The  variables  are  generally  treated  as 
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continuous  quantities.  If  discrete  variables  are  encountered,  they 
may  be  treated  initially  as  continuous  and,  upon  finding  an  optimum, 
a local  exploration  is  made  to  find  the  optimal  discrete  value  [12]. 


1.2  General  Strategy  of  Mathematical  Programming 

Like  all  optimization  methods  the  MP  methods  try  to  find  those 
values  x for  which  the  merit  function  f(x^)  will  be  minimized  (or 
maximized).  The  problem  may  be  stated  as  follows: 


Find  those  x^  that  produce  the 
min.f(x.. ) 

such  that  all  the  constraints 
9j(x.)  * 0 

and  the  equality  constraints 


j = 1 , 2 , . . . J 


(1.7) 


(1.8) 


0.9) 


gk(xi)  =0  k = J + 1,  J + 2,...K  (1.9) 

are  satisfied,  where  the  x^  are  the  I real  valued  variables  x-| , 
X2»...Xj  in  I-dimentional  Euclidian  space,  and  constraints  g^, 
g2,...gK  are  real-valued  real  functions  in  that  space. 

Most  MP  methods  do  not  treat  the  equality  constraints  directly. 
These  constraints  may,  for  example,  be  converted  to  inequality  con- 
straints of  the  form 

e ^ 9k(x1 ) ^ e (1.10) 

where  e is  an  arbitrary  small  number.  One  can  also  treat  an  equality 
constraint  by  solving  for  one  variable  in  terms  of  the  others,  thereby 
eliminating  a variable  and  constraint. 
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MP  methods  are  based  on  searching  strategies.  The  search  usually 
starts  from  an  arbitrarily  selected  point  x^  and,  by  means  of  some 
local  search  strategy,  a set  of  points  xr  is  generated  such  that 
f(xr)  < f(xr_1)  while  also  satisfying  all  the  constraint  equations. 

In  the  majority  of  design  problems,  either  f and/or  some  or  all  of 
the  constraints  g.  are  nonlinear  in  x,  and  therefore  one  has  a non- 
linear,  constrained  optimization  problem. 

None  of  the  available  nonlinear  optimization  methods  guarantees 
an  optimal  solution  (when  it  exists).  There  are  two  major  difficul- 
ties. First,  many  of  the  design  problems  are  not  unimodal  that  is, 
they  have  more  than  one  local  optimum.  The  nonlinear  methods  are, 
however,  designed  only  to  locate  local  optima  and  thus  the  global 
optimum  may  not  be  attained.  Secondly  the  search  algorithm  may  simply 
fail  to  find  even  a local  optimum. 

The  prudent  designer  thus  usually  tries  several  different  star- 
ting points.  If  all  the  search  paths  terminate  at  the  same  point, 
then  optimality  is  assumed  and  the  problem  can  be  considered  to  be 
unimodal.  However,  if  the  termination  point  is  different  for 
different  search  paths,  then  either  the  best  design  is 
accepted  or  additional  starting  points  are  tried  in  an  attempt  to 
find  a still  better  design.  Such  a decision  draws  on  the  designer's 
experience  and  judgment.  The  failure  to  converge  to  the  same  point 
may  be  due  either  to  the  presence  of  several  local  optima,  or  to 
algorithm  failure.  Failure  is  said  to  occur  when  the  algorithm  ter- 
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minates  at  a point  significantly  away  from  the  nearest  local  optimum. 
An  excellent  discussion  of  the  difficulties  involved  in  the  solution 
of  unconstrained  nonlinear  programming  problems  is  provided  in 
Reference  [12].  These  difficulties  are  greatly  compounded  when  con- 
straints are  added.  For  convenience,  since  this  and  other  MP  methods 
are  capable  only  of  locating  a local  optimum,  in  the  discussion  below 
the  term  optimum  should  be  taken  to  mean  local  rather  than  global 
optimum. 

1.3  Motivation  and  Structure 

Many  efficient  algorithms  exist  for  treating  linear  problems, 
that  is,  problems  with  linear  objective  functions  and  constraint 
equations,  and  certain  simple  nonlinear  problems  [5-7,  13].  Most 
design  problems  are,  however,  nonlinear  problems  that  must  be  treated 
by  relatively  less  efficient  methods. 

| 

There  are  a number  of  methods  available  for  solution  of  such 
problems  [14].  However,  there  are  several  difficulties  with  these 
methods.  For  example,  the  computations  required  to  produce  improved 
designs  in  the  neighborhood  of  the  constraint  boundaries  are  in  some 
cases  lengthy  and  relatively  complex.  The  number  of  function  evalua- 
tions is  often  very  large,  and,  most  methods  are  not  reliable,  that 
is,  the  best  points  produced  by  these  methods  may  not  be  a local 
optima. 


Recently,  Eason  and  Fenton  [15]  compared  a group  of  seventeen 
numerical  optimization  methods.  This  group  of  algorithms  contains 
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most  of  the  currently  popular  procedures.  The  results  of  this  study 
indicate  that  the  direct  search  algorithms  are  superior  with  respect 
to  generality,  efficiency,  running  cost,  speed,  preparation  cost  and 
reliability  [15].  Unfortunately,  none  of  the  methods  presented  in 
this  study  are  totally  reliable.  None  of  the  codes  solved  all  of 
the  ten  problems  attempted. 

In  a recent  paper,  Pappas  and  Moradi  compared  a code  based  on 
the  Direct  Search-Feasible  Direction  Algorithm  (DSFD)  [16]  with  those 
studied  by  Eason  and  Fenton  [17].  This  code  solved  all  the  test 
problems  treated  in  Reference  [15].  In  addition  it  also  solved  a 
difficult  six-variable  problem  which  the  relatively  reliable  DSDA, 
and  popular  SUMT,  procedures  failed  to  solve  [17].  In  this  study  the 
DSFD  based  code  proved  superior  to  all  others  in  overall  generality 
and  efficiency. 

Even  the  relatively  efficient  DSFD  however,  required  several 
hundred  to  several  thousand  function  evaluations  to  achieve  a solution 
to  Eason's  and  Fenton's  test  problems.  Therefore  substantial  compu- 
tational effort  and  cost  would  be  required  by  these  procedures  on 
problems  with  computationally  demanding  functions.  Since  many 
important  engineering  problems  are  of  this  type,  there  is  a clear 
need  for  a new  algorithm  which  requires  a reduced  number  of  function 
evaluations  for  solution. 


This  thesis  presents  an  apparently  fast  and  reliable  algorithm 
which  utilizes  a search  strategy  designed  to  greatly  reduce  the  number 
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of  function  evaluations  required  for  solution.  This  is  accomplished 
by  developing  a scheme  which  keeps  the  search  confined  to  the  neigh- 
borhood of  the  behavior  boundary,  on  the  premise  that  this  is  where 
the  optimum  is  usually  found  in  constrained  problems. 

The  search  moves  efficiently  along  this  boundary  by  making  effec- 
tive use  of  a relatively  small  amount  of  new  local  function  information 
relying  primarily  on  information  generated  earlier  in  the  search. 

Thus,  the  number  of  new  function  evaluations  required  to  sustain 
movement  and  thereby  the  total  number  of  function  evaluations  required 
for  solution  are  kept  low. 

Since  the  advent  of  high  speed  digital  computers  many  accurate 
but  computationally  demanding  methods  have  been  developed  for  engi- 
neering analysis.  Due  to  the  relatively  large  execution  time  required 
for  each  reanalysis  cycle  (function  evaluation)  of  many  of  these 
techniques  and  the  large  number  of  reanalysis  cycles  required  by  most 
MP  procedures,  the  combination  of  such  analytic  methods  and  MP  pro- 
cedures may  be  very  costly  to  utilize  and  are  therefore  often  imprac- 
tical. This  difficulty  can  seemingly  be  minimized  by  use  of  the 
Boundary  Tracking  (BT)  algorithm  due  to  its  apparent  efficiency  in 
reducing  the  number  of  function  evaluations  and  therefore  total  execu- 
tion time  required  for  solution. 

Where  the  BT  algorithm  is  not  available  in  a compiled  form,  a 
simpler  algorithm  may  be  preferable  for  use  on  small  simple  problems. 

In  such  circumstances  the  extra  compilation  time  required  due  to  the 


greater  complexity  of  the  BT  algorithm  may  overshadow  the  time  saved 
by  a reduction  in  the  number  of  function  evaluations  required  for 
solution. 


A 
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CHAPTER  2 

DESCRIPTION  OF  GENERAL  STRATEGY  OF 
THE  BOUNDARY  TRACKING  ALGORITHM 

2.1  Introduction 

The  most  common  strategy  for  constrained  optimization  techniques 
is  first  to  move  to  the  IF  boundary  from  a starting  point  and  then  to 
move  along  this  boundary  to  the  optimum.  A typical  objective  function 
surface  is  illustrated  in  Figure  1. 

The  general  strategy  of  the  Boundary  Tracking  (BT)  algorithm  is 
to  first  locate  a point  on  the  behavior  boundary.  Once  the  behavior 
boundary  is  located  the  best  direction  for  movement  along  this  boun- 
dary is  determined.  Movement  along  this  boundary  is  then  continued 
until  the  optimum  point  is  obtained. 

Figure  2 presents  a general  block  diagram  of  the  BT  algorithm. 

2.2  Location  of  the  Behavior  Boundary 

The  search  starts  from  an  arbitrary  starting  point  x9.  If  the 
problem  has  behavior  constraints,  these  constraints  are  evaluated  at 
this  point.  If  the  starting  point  proves  to  be  in  the  feasible  region, 
it  is  designated  as  a base  point  (point  1 Figure  1).  A step  is  then 
taken  in  the  direction  Of  the  gradient  of  the  objective  function.  All 
the  constraints  are  evaluated  after  the  step.  If  the  new  point  is  in 
the  feasible  region  the  objective  function  is  evaluated  and  compared 
to  the  last  base  point.  If  the  function  has  improved, this  new  point 
is  called  a base  point  (points  2,  3,  4 in  Figure  1)  and  the  search 


continues,  doubling  the  step  size  with  each  move. 


Once  the  behavior  boundary  is  crossed,  a point  on  this  boundary 
is  located  (point  5)  using  a line  search  in  the  direction  of  the  last 
move. 


The  search  is  performed  in  the  direction  of  movement,  rather 
than  say  the  direction  of  the  gradient  of  the  constraint  function 
violated,  since  identification  of  the  latter  direction  requires  cal- 
culation of  constraint  function  derivatives,  while  the  former  direc- 
tion is  available  and  thus  requires  no  new  information.  The  line 
search  is  used  because  of  its  simplicity. 

If  the  search  move  crosses  the  IF  boundary  at  a regional  limit, 
it  is  possible  to  find  the  best  direction  to  continue  movement  by 
solving  Zoutendijk's  direction  finding  problem  (see  section  3.2.3). 
Here  for  simplicity,  however,  the  variables  exceeding  the  limits  are 
set  equal  to  the  corresponding  limit  and  the  search  continued  until  a 
behavior  constraint  is  violated.  This  strategy  tends  to  deflect  the 
direction  of  movement  along  the  IF  boundary. 

If  after  a step,  the  point  remains  in  the  feasible  region  but 
the  objective  function  has  ^ot  improved,  a new  gradient  direction  is 
calculated  at  the  last  base  point  and  the  search  restarted.  If 
several  direction  changes  fail  to  locate  the  constraint  boundary,  a 
more  efficient,  albeit  more  complex,  unconstrained  search  method  such 
as  the  modified  Rotating  Coordinate  Pattern  (RCP)  search  [17],  is 


used.  This  search  procedure  is  applied  here  to  avoid  the  zig-zagging 
problem  associated  with  the  gradient  search  method  [12].  Thus  the 
advantages  of  the  gradient  search  method  are  exploited  where  possible, 
but  the  method  is  replaced  by  a more  appropriate  procedure  whenever 
strategy  dictates.  If  the  behavior  boundary  is  crossed  in  the  modi- 
fied RCP  search,  the  nearby  boundary  may  be  located  by  any  convenient 
method  such  as  that  already  described. 

The  major  contribution  of  this  thesis  is  the  procedure  for  move- 
ment along  the  behavior  boundary.  It  is  this  movement  which  causes 
the  difficulty  associated  with  the  solution  of  constrained  nonlinear 
problems.  The  method  used  for  locating  the  boundary,  such  as  the 
gradient  search  or  the  RCP  search,  is  of  secondary  importance.  Since 
the  bulk  of  the  search  is  associated  with  movement  along  the  boundary, 
the  choice  of  the  initial  boundary  locating  procedure  will  not  sub- 
stantially effect  overall  performance.  Thus,  the  procedures  used  for 
this  purpose  are  not  of  great  importance  and  need  not  be  described  in 
detail  here.  The  interested  reader  is  refered  to  references  [12,  17] 
for  a detailed  description  of  these  methods. 

If  the  starting  point  is  infeasible,  the  behavior  boundary  is 
located  by  first  selecting  the  most  negative  constraint,  that  is,  the 
constraint  presenting  the  greatest  violation.  The  point  where  this 
constraint  vanishes  is  then  located  using  a line  search  in  the  direc- 
tion of  the  gradient  of  this  constraint  (see  section  3.1.2).  This 
direction  is  used  here  since  it  provides  the  shortest  distance  to  the 
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surface  where  this  constraint  vanishes.  The  line  search  procedure 
used  here  is  a modified  Secant  root  finding  procedure  (see  section 
3.1.2). 

All  the  constraints  are  then  evaluated  at  the  point  where  the 
selected  constraint  vanishes.  If  no  violation  exists  the  objective 
function  is  evaluated  and  the  point  is  designated  as  a base  point. 
However,  if  any  other  constraints  are  violated  at  this  point,  the 
smallest  one  (algebraically)  is  again  selected  and  the  point  where  it 
is  zero  is  found.  This  process  is  continued  until  a point  on  the 
behavior  boundary  is  located. 

For  problems  without  behavior  constraints  an  unconstrained  search 
method,  such  as  the  modified  RC  Pattern  search  [18],  is  used  to  locate 
the  optimum  point.  This  search  method  combines  the  well  known  RC 
Pattern  search  with  Zoutendijk's  feasible  direction  finding  method 
(see  section  3.2.3).  The  direction  finding  procedure  is  employed  at 
the  points  of  RC  Pattern  search  failure. 

2.3  Movement  Along  The  Behavior  Boundary 

Once  the  point  on  the  behavior  boundary  is  located  it  is  designa- 
ted as  a base  point  (point  in  Figure  1).  The  best  direction  to 
move  is  then  determined  by  solving  the  direction  finding  problem  of 
Zoutendijk  (see  section  3.2.3)  and  a step  taken  in  this  direction. 

The  objective  function  is  then  evaluated  at  this  point.  If  the  func- 
tion has  improved,  all  the  constraints  are  evaluated. 
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If  none  of  the  constraints  that  were  inactive  at  the  base  are 

violated,  an  appropriate  constraint  is  selected  from  among  the  active 

constraints  (see  section  3.1.1).  The  behavior  boundary  is  located 

by  finding  a point  where  this  constraint  is  zero  by  means  of  a line 

search  in  the  direction  of  the  gradient  of  the  selected  constraint 
o 

(point  x ). 

Upon  locating  the  point  on  the  boundary  after  the  best  direction 
move,  the  objective  function  is  evaluated.  If  the  function  has  im- 
proved  the  point  is  designated  as  a base  point  (point  x ).  A move  is 
then  made  from  the  last  base  in  the  direction  of  improvement  along  a 
line  connecting  two  bases  a distance  equal  to  the  distance  between 
these  bases. 

The  boundary  is  then  located  as  before  (see  point  x of  Figure  1) 
along  the  direction  of  the  gradient  of  the  appropriate  constraint  as 
evaluated  at  the  last  point  where  the  best  direction  finding  problem 
was  formulated.  It  should  be  pointed  out  that  this  is  an  underlying 
feature  of  this  algorithm,  making  multiple  use  of  function  evaluations 
in  order  to  keep  the  number  of  new  evaluations  required  as  small  as 
possible.  Only  if  this  direction  fails  to  locate  the  boundary,  are 
new  gradient  values  computed. 

When  the  point  on  the  behavior  boundary  is  found,  it  is  compared 
to  the  last  base  point.  If  the  function  has  improved,  the  point  is 
designated  as  a new  base  point  (points  x , x , x ) and  the  search  is 
continued.  However,  if  the  function  has  not  improved  the  direction 


1 


is  abandoned  and  the  search  is  re-started  from  the  last  base  point  by 
calculating  a new  best  direction. 

After  each  move  all  the  constraints  are  evaluated,  if  a new  con- 
straint is  violated  (point  A)  the  point  where  this  constraint  is  zero 
is  located  and  a new  best  direction  calculated.  However,  if  no  new 
constraint  is  violated,  an  appropriate  constraint  is  selected  from 
among  the  previously  active  constraints  (see  section  3.1.1)  and  the 
boundary  located,  the  process  is  continued  until  convergence  is 
achieved. 

If  no  appropriate  constraint  is  active  the  move  along  the  pre- 
vious direction  is  continued,  since  it  is  the  best  direction  avail- 
able. If  after  a step  in  the  above  direction,  the  objective  function 
has  increased  rather  than  decreased  the  step  size  is  reduced.  The 
process  is  repeated  until  a better  point  is  found  or  the  convergence 
criterion  is  met  (see  section  3.3). 

If  after  any  move  a constraint  inactive  at  the  last  base  point 
is  violated,  the  behavior  boundary  is  located  by  finding  the  point 
where  this  constraint  is  zero.  This  point  is  found  by  a line  search 
in  the  direction  of  the  last  move.  This  direction  is  used  to  elimi- 
nate the  need  to  calculate  the  gradient  of  the  new  constraint.  The 
search  is  then  re-started  by  calculating  the  best  movement  direction 
at  this  point  on  the  boundary. 


2.4  Search  termination  description.  For  the  few  cases  where 
the  optimum  occurs  away  from  the  behavior  boundary,  a sufficient  con- 
dition for  optimality  is  that  all  components  of  the  gradient  of  the 
objective  function  be  essentially  equal  to  zero.  For  most  cases, 
however,  where  the  optimum  is  on  the  behavior  boundary  no  such  simple 
condition  exists. 

The  search  termination  procedure  used  here  is  as  follows.  The 
search  using  a specified  basic  step  size  is  performed  until  no  further 
improvement  can  be  obtained.  When  such  a point  is  found  it  is  desig- 
nated as  an  optimality  comparison  base.  Now  the  basic  step  size  is 
reduced  in  an  effort  to  achieve  improvement  and  the  search  continued 
until  a new  optimality  comparison  base  is  obtained. 

In  order  to  justify  a further  reduction  in  the  step  size  an 
optimality  check  is  performed  whenever  the  step  size  is  to  be  reduced. 
If  a convergence  limit  is  satisfied,  the  step  size  is  again  reduced 
for  further  improvement  and  the  search  is  continued  until  another 
optimality  comparison  base  is  obtained.  A secondary  convergence 
check  is  then  performed  in  order  to  determine  whether  the  latest  im- 
provement due  to  the  step  reduction  is  less  than  the  previous  improve- 
ment thereby  indicating  convergence  has  occurred.  If  both  convergence 
tests  are  simultaneously  satisfied,  or  the  minimum  step  size  is  reach- 
ed, the  search  is  terminated.  Otherwise,  the  search  is  re-started  by 
calculating  a new  best  movement  direction. 
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A block  diagram  of  the  search  termination  procedure  is  given  in 
Figure  3. 


Figure  3 - Block  diagram  of  the  termination  procedure 
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MATHEMATICAL  DEVELOPMENT 
3.1  Boundary  Location 

It  is  assumed  for  the  purpose  of  devising  a search  strategy  that 
the  solution  to  the  optimization  problem  with  behavior  constraints 
lies  on  the  behavior  boundary.  The  main  strategy  is  first  to  locate 
such  a boundary  and  then  to  move  along  this  boundary.  For  the  pur- 
pose of  locating  the  boundary  a "proper"  behavior  constraint  is  iden- 
tified and  method  for  locating  the  behavior  boundary  selected  or 
developed. 


3.1.1  Identification  of  "proper"  constraints.  All  negative  and 
small  positive  constraints  are  designated  as  active  constraints  (for 
definition  of  activity  see  section  3.2.3).  A constraint  is  considered 
"proper"  if  it  is  necessary  or  desirable  to  locate  a point  where  the 
value  of  this  constraint  is  zero  after  making  a search  move  in  an 
effort  to  establish  a new  base. 


Since  it  is  necessary  to  eliminate  any  constraint  violation  pre- 
sent at  a point,  all  the  negative  constraints  are  considered  "proper". 

If  all  active  constraints  are  positive,  a positive  active  con- 
straint is  considered  proper  if 

- Vf  • Vg.  < 0 (3.1) 

J 

This  indicates  that  the  projection  of  -Vf  on  Vg  is  of  opposite 
sign  to  Vg  (point  1 Figure  4)  and  therefore,  the  value  of  the  constraint 


- 


i 
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will  decrease  as  a result  of  a move  in  -Vf  direction  (minimization 
problem) . 

The  best  direction  for  improving  a function  is  its  gradient 
direction.  However,  this  direction  is  not  always  feasible  in  con- 
strained problems  (point  1 Figure  4). 

A direction  is  desired  such  that  after  taking  a step,  the  objec- 
tive function  will  improve  the  maximum  amount  possible  without 
violating  the  behavior  constraints.  A direction  along  the  behavior 
boundary  often  has  this  property.  Moving  along  the  behavior  boundary, 
however,  may  not  always  produce  the  best  improvement  in  the  objective 
function  as  shown  by  point  2 in  Figure  4.  Here  projection  of  -Vf  on 
Vg.  has  the  same  sign  as  Vg..  Thus,  the  value  of  the  constraint  g. 

J J J 

will  increase  by  moving  in  -Vf  direction.  A criteria  is  therefore 
needed  that  can  indicate  which  positive  constraint  (if  any)  is  "pro- 
per". Equation  (3.1)  satisfies  this  purpose,  since  if  satisfied  it 
identifies  the  constraints  that  may  be  violated  if  a move  in  the  best 
direction  (-Vf)  is  made.  Thus  it  identifies  proper  constraints, 
that  is,  those  constraints  along  which  it  is  desirable  to  move  in 
order  to  achieve  the  best  improvement  in  the  objective  function. 

If  no  proper  constraint  can  be  identified,  the  movement  along  the 
previous  direction  is  continued,  since  there  is  no  advantage  in  loca- 
ting and  moving  along  the  nearby  boundary. 

3.1.2  The  boundary  locating  method.  Call  the  set  of  proper 
constraints  Pk  where,  k = 1,  2,...K.  If  more  than  one  "proper" 


1 


Figure  4 - Selction  of  a Proper  Constraint 
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constraint  is  identified,  that  is  if 

K > 1 (3.2) 

the  smallest  (algebraically)  is  selected  and  called  g*.  The  root  of 
g*  is  then  found  in  the  direction  of  movement  or  the  gradient  of  g*. 

The  root,  if  it  exists,  of  the  constraint  g*(x)  along  a line  may 
be  found  as  closely  as  desired  by  means  of  the  "secant  root-finding" 
method  [19]  by  iterating  the  equation 

xs+1  = xs  - g*(xs)(xs-xs_1 )/[g*(xs)  - g*(xs_1)]  s = 1,2, (3.3) 
where  N1  is  the  maximum  number  of  iterations  permitted,  until 

|g*  (xS+1)|  (3.4) 

is  satisfied.  Here  e-j  is  an  arbitrarily  selected  accuracy  limit. 

S c — I 

Initially  x , x are  two  points  on  the  line  defined  by 
Case  A:  If  the  boundary  is  to  be  located  in  the  direction  of  the  last 

search  move  (see  section  2.3)  let  the  initial  points  on  the  line 
(s  = 1)  be  given  by 

(3.5) 

(3.6) 


x5-1  - xb 


xs  = xr 


b r 

where  x and  x are  the  end  points  of  the  last  search  step. 

Case  B:  If  the  boundary  is  to  be  located  in  the  gradient  direction 

of  the  constraint  g*(x)  (see  section  2.3)  let 

s-1  r 
x = x 


xs  = xr  ± anVg*  (xn)  / (vg*  (xn) ) 
joint  from  which  the  boundary  is  tc 
is  the  last  base  where  a direction  finding  problem  was  formulated 


(3.7) 

(3.8) 


where  xr  is  the  point  from  which  the  boundary  is  to  be  located  and  xn 


/ 
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(see  section  3.2.3.).  Here  the  plus  and  minus  signs  are  used  when  g* 
is  negative  or  positive  respectively.  The  quantity  an  is  the  step  size 
used  at  point  x11. 

Case  C:  For  the  case  where  this  procedure  is  used  to  find  the  behavior 

boundary  where  a starting  point  of  the  optimal  search  is  infeasible, 

v*  Q 

replace  x by  the  starting  point  x and  let, 

a"  = a°  (3.9) 

where  a0  is  an  arbitrarily  specified  initial  step  size,  and  let 

Vg*  (xn)  = Vg*  (x°)  (3.10) 

in  equations  (3.7)  and  (3.8). 

After  equation  (3.4)  is  satisfied,  the  remainder  of  constraints 
in  the  set  may  also  vanish.  However,  if  after  eliminating  the 
the  smallest  violation 

gj(xs+1)*-e1  j = 1 ,2 ...  .J  (3.11) 

is  not  satisfied,  a new  g*  is  selected  and  the  procedure  repeated  until 
a point  on  the  behavior  boundary  is  located  (equation  3.11  satisfied). 

In  case  B,  slow  convergence  may  indicate  that  the  gradient  Vg*(xn) 
may  no  longer  be  applicable  at  point  xn.  This  situation  is  shown  in 
Figure  (5).  There  is  no  point  along  the  gradient  Vg*(x^)  which  falls 
on  the  boundary.  A new  gradient  direction  Vg*(x  ) is  therefore  calcu- 
lated and  used  in  place  of  Vg*(xn).  If  after  calculating  a new  gradi- 
ent direction  the  boundary  still  cannot  be  located  with  a reasonable 


number  of  iterations,  the  difficulty  may  be  due  to  too  large  a reach 
6 r 

step  (line  x - x in  Figure  5).  Thus  where  recalculation  of  the 
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gradient  fails,  the  step  size  is  reduced  and  another  effort  is  made  to 
locate  the  boundary  from  a point  somewhat  closer  to  the  last  base. 

Therefore,  unless 

s < N1  . (3.12) 

where  N1  is  chosen  empirically,  let 

Vg*  (xn)  = vg*  (xr)  (3.13) 

in  equation  (3.8)  and  repeat  the  application  of  equation  (3.3)  If  equa- 
tion (3.12)  is  again  violated  let 

Ax1'*1  = Axr/2  (3.14) 


where 


r r b 
Ax  = x - x 


(3.15) 


and  x is  the  last  base  point.  Then  select  the  smallest  proper  con- 

P+1 

straint  at  point  x where 

x^1  = xb  + Axr+1  (3.16) 

f 1 y,  M 1 1 

(point  x + Ax  in  Figure  5),  replace  x by  x in  equations  (3.7) 
and  (3.8)  and  repeat  the  iteration  of  equation  (3.3).  The  process  is 
continued  until  the  boundary  is  located  without  violating  equation 
(3.12)  or  until 


|Axr+1|  <an.  (3.17) 

If  equation  (3.17)  is  satisfied  the  attempt  to  locate  the  boundary  is 
abandoned.  The  optimal  search  is  then  restarted  by  calculating  a new 
starting  direction  (see  section  3.2.3)  at  the  last  feasible  point  en- 
countered. 


3.2  Mathematical  Representation  of  General  Strategy 

The  general  strategy  of  the  algorithm  can  be  represented  mathe- 
matically as  follows. 


28. 


3.2.1  Movement  to  the  boundary-starting  point  in  the  feasible 
regi on . A block  diagram  for  the  movement  to  the  boundary  procedure 
is  given  in  Figure  6. 


Given  the  starting  point  xC,  where  initially  r = 0,  which  satis- 


fies regional  constraints 


and  all 


L r U 
x.  ;<  x.  ^ xi 


9j<x1>  *° 


i = 1,  2,... I 


j = 1,  2....J 


(3.18) 


(3.19) 


the  objective  function  f(xr)  is  evaluated.  This  point  is  called  a 
base  point  xb,  where  initially  b = n = 1 , and  the  first  comparison  base 


xn  is  defined  as. 


xn  = xb 


(3.20) 


Then  with 


r _ n 
a = a 


(3.21) 


xr+l  = xr  - ar  [Vf ( xn ) ] / | Vf (x11) 


(3.22) 


I 

where  V is  the  gradient  operator,  | Vf | the  magnitude  of  Vf  and,  a is 
the  step  size  at  iteration  r. 


If  any 


r+1  U 
xi  * xi 


r+1  U 
xi  = xi 


(3.23) 


or  if  any 


J 


Xr+1  <XL 

i xi 


r+1  L 
xi  s xi 


(3.24) 


Now  compute  the  g^(x  ),  if  equations  (3.19)  are  satisfied,  evaluate 
f(xr+1).  If 


f(xr+1)  < f(xb) 

call  the  point  a new  base  point.  Thus  let 

b+1  r+1 

x = x 


(3.25) 


(3.26) 


A new  move  is  then  made  according  to  equation  (3.22)  with  the 


step  size  redefined  as 


r+1  „ r 

a = 2a 


(3.27) 


and  application  of  equations  (3.23-3.26)  repeated.  If  b>l  and  equa- 
tion (3.25)  is  not  satisfied  increase  n by  one,  let  the  next  compar- 
ison base  xn  = xb  (3.28) 


ison  base 
and  let 


Now  if 


r b 
x = x 


n < N, 


(3.29) 


(3.30) 


where  N2  is  an  arbitrary  constant  signifying  the  number  of  direction 
changes,  a new  gradient  direction  Vf  is  calculated.  The  application 
of  equations  (3.20-3.27)  are  then  repeated  until  either  a point  in 
the  infeasible  region  is  found  or  equation  (3.30)  is  not  satisfied. 

r+l 

If  point  x is  in  the  infeasible  region  the  behavior  boundary  is 
located  in  the  direction  of  the  last  move  (see  section  3.1.2). 


If  n exceeds  N2  the  gradient  search  is  abandoned  and  the  RCP 
search  invoked.  This  search  method  is  continued  in  the  unconstrained 
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region  until  a point  in  the  infeasible  region  is  found  or  the  optimum 
is  achieved  (see  section  3.3). 

If  equation  (3.25)  is  not  satisfied  when  b=l , the  step  size  a i;s 


redi fined  as 


r+1  r/0 

a = a /2 


(3.31) 


and  a new  move  is  made  according  to  equation  (3.22).  The  application 
of  equations  (3.31)  and (3.22)  is  repeated  until  equation  (3.25)  is 


satisfied  or 


r+1  ^ 

a < “min 


(3.32) 


where  is  a minimum  step  size.  In  the  latter  case  the  point  is 
assumed  to  be  optimum  and  the  search  is  terminated. 


3.2.2  Movement  to  the  boundary  - starting  point  in  the  infeas- 
ible region.  If  the  starting  point  x°  is  in  the  infeasible  region  or 
the  RCP  search  has  located  a point  in  this  region  the  behavior  bound- 
ary is  located  by  a search  in  a constraint  gradient  direction  (see 
section  3.1.2). 

3.2.3  Initiation  of  movement  along  the  boundary.  After  locating 
a point  on  the  behavior  boundary  a method  is  needed  which  can  initiate 
the  move  along  this  boundary  toward  the  optimum  if  optimality  has  not 
yet  been  achieved.  The  Feasible  direction  finding  algorithm  of 
Zoutendijk  [20]  is  suitable  for  this  purpose  since,  it  either  provides 
a direction  along  which  an  improved  point  can  be  found,  or  indicates 
the  presence  of  a local  optimum. 


Definition:  a)  A direction  vector  u is  called  "feasible"  if 

after  taking  a sufficiently  small  step  along  this  direction  no  con- 
straint is  violated  [20].  This  will  be  true  if 

(u)T  Vg^  (x)  » 0 (3.33) 

where  (u)T  is  the  transpose  of  u,  since  a small  step  in  this  direc- 
tion will  produce  no  change  or  an  increase  in  g.(x). 

J 


b)  A feasible  vector  u is  called  "usable"  if  a move  in.  the  direc- 
tion u can  also  provide  an  improvement  in  f(x).  This  will  be  the  case 
if 

(u)T  Vf (x)  < 0 (3.34) 


since  a sufficiently  small  step  in  the  u direction  will  produce  a 
decrease  in  the  value  of  f(x). 

c)  A behavior  constraint  g.(x)  is  considered  "active"  if 

J 

9j(xi ) « (3.35) 

where 

e.  = K,  (an  / a0)  (3.36) 

J J 

is  an  array  of  positive,  arbitrary  small  numbers.  These  numbers  may 
approach  zero  during  the  search,  but  can  not  be  identically  zero  [20]. 
The  quantity  a*  is  the  step  size  at  the  beginning  of  the  search  and 
K.  is  a small  positive  constant. 

J 

d)  A lower  regional  constraint  is  considered  "active"  if 

x1-x![4a n (3.37) 

and  the  upper  constraints  active  if 

x“  - x.  <:  a"  (3.38) 


| 


33. 


Denote  the  set  of  all  active  behavior  constraints  by  A..  Call 

J 

R.j  and  the  active  lower  and  upper  regional  constraint  set  respec- 
tively. 

The  direction  finding  problem  can  be  formulated  in  the  following 
manner: 


Given  x,  find  u that  results  in  the 


max  a 

(3.39) 

and  for  which 

a > 0 

(3.40) 

(u  )T  Vf  (x)  + a ^ 0 

(3.41) 

-(u  )T  Vgj(x) 

+ W-a^0  J € A . 

J J 

(3.42) 

ui  ^ 0 

i € Rj 

(3.43) 

u.j  » 0 

i € R\ 

(3.44) 

M $ 1 

i = 1 , 2 ....  I 

(3.45) 

Equations  (3.45)  bound  the  length  of  u to  prevent  the  direction 
finding  problem  from  producing  an  unbounded  solution  vector. 

When  suitable  deflection  parameters  W.  are  used,  the  solution  of 

J 

the  problem  provides  a usable  and  feasible  direction  u since,  if  a > o 
from  equations  (3.41)  and  (3.42)  equations  (3.33)  and  (3.34)  will  be 
satisfied.  It  also  provides  the  best  usable  direction  since,  if  o is 
maximized,  the  left  hand  side  of  equation  (3.34)  will  be  a maximum 
and  therefore  the  direction  found  will  be  the  best  direction  for  de- 
creasing f(x)  [20]. 
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It  may  be  seen  that  the  direction  finding  problem  of  equations 
(3.39  - 3.45)  is  a linear  programming  problem.  Thus,  any  linear 
programming  method  such  as  the  simplex  procedure  [21]  can  be  used  to 
provide  the  solution. 

Since  the  objective  here  is  to  move  along  the  IF  boundary,  the 
direction  u should  be  as  close  to  this  boundary  as  possible.  The 
work  associated  with  locating  the  boundary  will  therefore  be  minimized 
after  the  move.  The  deflection  parameters  W.  here  are  thus  set  to  be 

J 

equal  to  zero.  This  produces  a direction  u that  tends  to  be  tangent 
to  the  boundary. 

3.2.4  Movement  along  the  boundary.  A flow  chart  for  the  move- 
ment along  the  boundary  procedure  is  given  in  Figure  7. 


Call  the  point  on  the  behavior  boundary  xD  and  define 


a1  = an 


(3.46) 


where  initially  b = l - 1,  xb  is  a base  point,  and  c£is  a comparison 
step  size.  A direction  u is  then  calculated  (see  section  3.2.3)  and 


the  point  xr  defined  as 


where 


x"  = x*  ♦ A xj 


Ax'  = a*  u.  / j u | 


(3.47) 


(3.48) 


f(xr)  < f(xD) 


(3.49) 


9j(x  ) > e. 


j = 1,  2....J 

j *Aj 


(3.50) 


■ 35. 
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point  on  the 
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ncr.  I by  one. 
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direction,,  v. 
u,  let  f =f(x  ) 


Optimality 

Procedure 


/n  , 
a <a  . ? 
mm 


yI  t u/  Point  x 

n o r\ Optimum?. 


Identify  set  | [Let  aJl_  n 


'Aj  chgd . 


r b l . 
xi=xi+a  u^/  u 


vi<f(xV^  g 


Locate  the  ^ 
boundary  in  Ax 
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e.  (xr-)>e.?\  N Locate  the  /T* 

J ^ boundary  in  the 

v ^rAj‘  direction  u 
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direction 


T\X^CxS)<f(xb)? 


Figure  7 - Flow  Chart  of  the  Procedure  for  Movement  along  the  Boundary 


, 


are  satisfied,  point  xs  on  the  behavior  boundary  is  located  by  a 
search  in  the  constraint  gradient  direction  and  the  objective  func- 
tion evaluated.  Now,  if  equation  (3.49)  is  satisfied  at  xs,  let 

xb+1  = xs  (3.51) 


now  index  b and  define  the  next  step  as 


r b r 

Xi  = xi  + A xi 


where 


i Jj  = • x^2 


4 x(  * X1  - X?'1 


b > 2 


b = 2 


(3.52) 


(3.53) 


(3.54) 


The  constraints  are  now  evaluated.  If  equation  (3.50)  is  satisfied, 
the  behavior  boundary  is  located  in  the  constraint  gradient  direction 
and  the  application  of  equations  (3.51  - 3.54)  repeated. 

If  when  b=l,  equation  (3.50)  is  not  satisfied,  that  is,  if  new 
constraints  have  become  active,  the  value  activity  limit  e.  is  doubled 

J 

for  all  constraints  in  set  A.  not  satisfying  equation  (3.50),  since 

J 

its  previously  assigned  value  was  not  sufficient  to  properly  define 
activity.  The  boundary  is  then  located  in  the  direction  u and  the 
application  of  equations  (3.47  - 3.50)  repeated. 

If  equation  (3.50)  is  not  satisfied  when  b > 2 the  boundary  is 
located  in  the  direction  of  the  last  move  Axr.  Then  a new  set  of 
active  constraints  is  defined  and  the  search  re-started  from  the  last 
base  point  (see  section  3.2.3). 
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In  application  of  equations  (3.49  - 3.54)  if  equation  (3.50)  is 
satisfied,  but  equation  (3.49)  is  not  satisfied,  this  indicates  that 
the  last  search  move  direction  was  unproductive,  therefore  the  last 
direction  of  movement  is  abandoned  and  the  search  is  re-started  from 
the  last  base  point  xb  by  calculating  a new  direction  u,  and  setting 
b*l . 


After  employing  equations  (3.47)  and  (3.48),  if  equation  (3.49) 
is  not  satisfied,  this  indicates  that  the  step  size  is  too  large, 
therefore  let 


n+1 


an/2 


(3.55) 

(3.56) 

Now,  since  new  activity  limits  are  defined  [see  equation  (3.36)]  the 


£+1  _ n+1 
a = a 


active  constraints  are  again  determined,  if  the  set  A_-  has  changed  a 

J 

new  direction  u is  calculated,  otherwise  the  previous  direction  is 
used  and  application  of  equations  (3.47  - 3.49)  repeated  using  the 
new  step  size.  This  procedure  is  continued  until  equation  (3.49)  is 
satisfied  or  convergence  is  achieved. 


Application  of  equations  (3.47  - 3.56)  is  invoked  when  applicabl 
as  explained  above,  until  the  convergence  or  optimality  criterion  are 
met  (see  section  3.3). 


3.3  Search  Termination 

A flow  chart  for  the  search  termination  procedure  is  given  in 


Figure  8. 
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The  question  of  optimality  is  considered  at  all  points  where  the 
direction  finding  problem  is  formulated.  For  the  few  cases  where  a 
local  optimum  occurs  away  from  the  behavior  boundary  that  is  where 
all 


9j  (X)  > £j  (3.57) 

a sufficient  condition  for  optimality  is 

|Vf|  « e2  (3.58) 

where  e2  is  an  arbitrary  small  number. 

For  most  cases,  however,  where  optimum  is  constrained,  the  solu- 
tion to  Zoutendi jk's  direction  finding  problem  [20]  provides  a test 
for  optimality.  A null  solution  vector  u indicates  a local  optimum 
where  all  e.  = o.  However,  where  some  e.  + o and  the  active  con- 

J J 

straints  are  also  not  zero,  the  point  may  be  merely  near,  rather  than 
at  the  optimum,  since  in  this  type  of  problems  the  optimum  point 
is  usually  on  the  behavior  boundary.  Therefore,  when  u = o attempts 
are  made  to  re-start  the  search  and  the  optimality  procedure  is  invoked 
only  when  a non- zero  direction  vector  u is  obtained.  In  order  to  re- 
start the  search  the  step  size  is  redefined  as 

an+1  = an/2  (3.59) 

and  therefore  the  activity  limits  Cj  are  redefined  according  to  equa- 
tion (3.36).  A new  direction  u is  then  calculated  if  A.  has  changed. 

J 

This  procedure  is  repeated  until  a non-zero  direction  is  found  or  the 
minimum  step  size  is  encountered,  that  is 

^ °Win  (3-60> 

in  which  case  the  search  is  terminated  and  the  point  is  assumed  to  be 


i 
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an  optimum. 

For  the  cases  where  u f o,  increase  l by  one  and  let 

</  = an  (3.61) 

fA  = f(xb)  (3.62) 

o n 

where  f is  a comparison  value  with  initially  Z = 0,  and  f = C-|, 
where  C-|  is  an  arbitrary  large  number.  Now  if 

a1  = a*"1  (3.63) 

no  optimality  check  is  made  since  Z has  increased  because  the  search 

was  re-started  due  to  the  fact  that  the  previous  movement  direction 
was  unproductive.  However,  if  the  step  size  has  changed,  this  indi- 
cates that  the  maximum  improvement  in  the  function  has  been  achieved 
using  the  previous  step  size  and  the  step  size  must  be  reduced  for 
additional  improvement.  A convergence  check  is  now  made  in  order  to 
determine  if  a search  using  the  new  step  size  is  justified.  Thus, 
where  equation  (3.63)  is  not  satisfied,  define 

Af4  = |(f4  - f*"1)  / f*|.  (3.64) 

Unless 

Af4  < e3  (3.65) 

where  e^  is  the  primary  convergence  criteria,  the  search  is  continued. 

A secondary  convergence  check  is  initiated  whenever  the  primary 
convergence  criteria  is  met,  [equation  (3.65)  satisfied]  in  order  to 
confirm  optimality.  For  this  purpose  the  step  size  a is  reduced 
according  to  equation  (3.59)  and  the  search  is  re-started.  When 
further  movement  with  the  new  step  size  terminates  the  primary  conver- 
gence check  is  invoked.  Therefore,  whenever  equation  (3.65)  is  not 


L 


I 
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satisfied  at  point  xb,  let 

K = 0 

(3.66) 

and  where  equation  (3.65)  is  satisfied  index  K and  let 

an+1  = a"/2 

(3.67) 

£+1  _ n+1 

a = a 

(3.68) 

If  equation  (3.60)  is  satisfied,  the  point  is  assumed  to  be  a local 
optimum  and  the  search  is  terminated.  Otherwise  a new  direction  u is 
computed  and  the  search  continued.  When  the  primary  convergence  check 
is  satisfied  in  two  consecutive  tries,  that  is  if 

K = 2 (3.69) 

the  secondary  convergence  check  is  invoked.  Thus  if 

Af4  < Af4'1  (3.70) 

the  search  is  terminated.  Otherwise  a new  direction  u is  computed,  K 
is  set  equal  to  1 and  the  search  continued.  When  equation  (3.70)  is 
satisfied,  the  change  in  the  value  of  the  objective  function  (Af)  has 
decreased  from  the  previous  convergence  check,  even  though  a reduced 
step  size  was  used,  the  point  is  thus  assumed  to  be  optimum  and  the 
search  is  terminated. 
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CHAPTER  4 
DESIGN  EXAMPLE 

4.1  Problem  Statement 

Consider  the  minimum  weight  design  of  submersible,  circular, 
cylindrical  shells  reinforced  by  equally  spaced  "T"  type  frames.  The 
design  variables  with  respect  to  which  the  optimization  is  to  be  car- 
ried out  are:  Plate  thickness,  frame  web  and  frame  flange  thickness, 

frame  flange  width,  web  height,  and  frame  spacing.  The  objective 
function  to  be  minimized  is  the  ratio  of  shell  weight  to  the  weight 
of  fluid  displaced.  All  the  standard  design  equations  used  in  submers- 
ible shell  design  practice  are  to  be  satisfied.  All  variables  will 
be  treated  as  continuous.  The  fixed  design  parameters  are:  The 

operating  depth,  shell  diameter,  shell  segment  length,  shell  eccen- 
tricity, construction  material  properties,  factors  of  safety  to  be 
used  in  design,  maximum  and  minimum  values  permitted  for  the  design 
variables,  and,  when  required,  a maximum  (when  external  frames  are 
used)  or  minimum  (when  internal  frames  are  used)  frame  diameter. 

This  problem  is  selected  here  since  it  is  a difficult  and  compu- 
tationally demanding  engineering  problem  which  the  relatively  reliable 
optimization  code,  DSDA  and  the  popular  SUMT  procedure,  could  not 
solve  [17]. 

In  a previous  study  [22]  it  has  been  shown  that  the  merit  surface 


in  this  problem  is  fairly  flat,  therefore  a wide  range  of  variables 
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with  similar  objective  function  values  may  be  generated  during  the 
search.  This  study  also  demonstrates  that  movement  along  the  behavior 
boundary  is  very  difficult  in  this  problem. 


A modified  version  of  this  problem  with  only  four  variables  was 
treated  in  Reference  [22],  using  the  DSDA  procedure.  The  formulation 
of  the  more  difficult  six  variable  problem  presented  here  is  similar 
to  that  used  in  Reference  [22]. 

Only  CAD0P3  [23],  a modified  version  of  CAD0P2  code,  [16]  reached 
an  optimal  solution  to  this  problem.  This  code  is  three  times  faster 
than  CAD0P2  code  in  problems  with  behavior  constraints.  However,  due 
to  the  difficulties  involved  in  moving  along  the  behavior  boundary,  the 
execution  time  required  by  CAD0P3  to  achieve  the  optimum  is  quite  long 
(see  section  4.3).  The  problem,  therefore,  was  treated  using  the  BT 
algorithm  in  order  to  demonstrate  its  superior  ability  in  moving  along 
a difficult  behavior  boundary. 


4.2  Problem  Formulation 


4.2.1  Objective  function.  For  cylindrical  vessels  with  periodic 
"T"  type  reinforcing  frames  the  objective  function  may  be  written  as 


f(x)  = 


W^w  VD 


internal  frames 


W/[yw(VD  + Vw  + Vf)]  external  frames 


where  W is  the  weight  of  the  hull  segment,  given  by 

W - Y$  (V$  ♦ Vw  + Vf) 


(4.1) 


(4.2) 
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The  y_  and  v are  the  specific  weight  of  the  material  and  immersion 
s 'W 

fluid,  respectively.  Vs>  Vw>  Vf  are  the  volume  of  the  plating,  frame 
webs,  and  frame  flanges,  respectively,  and  VQ  is  the  displacement 
volume  of  the  cylinder  enclosed  by  the  plating  envelope.  The  varia- 
bles , Xg.-.-jXg,  which  in  turn  represent  the  quantities  t,  b,  b', 
w,  L1 , and  h,  respectively,  are  defined  in  Figure  9.  The  objective 
function  f(x),  therefore  represents  the  ratio  of  the  shell  weight  to 
the  weight  of  the  fluid  displaced  for  the  shell  segment  of  length  L$. 
The  weight  of  the  bulkheads  is  omitted. 

4.2.2  Constraint  equations.  The  constraint  equations  g.(x)  are 

J 

divided  into  two  groups:  a)  Behavior  constraints  which  control  the 

failure  modes  or  impose  limitations  on  the  space  relationships  among 
the  variables,  and  b)  regional  constraints  which  specify  ranges  of 
the  variables.  The  basic  behavior  constraint  equations  are  formulat- 
ed using  a modified  version  of  the  design  equations  given  in  References 
[24]  ans  [25]. 

The  general  instability  constraint  is  given  by 

gl  = (P*cg  - S-,p)  / P*cg  * 0 (4.3) 

where  p is  the  applied  hydrostatic  pressure;  p*cg  is  the  minimum  of 

the  collapse  pressure  p (n,m)  due  to  general  instability  [26],  and 

cy 

S-j  is  the  factor  of  safety  for  this  failure  mode.  The  collapse  pres- 
sure is  given  by 


Figure  9 - Typical  Shell  Cross  Section 
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Pcg(n,m)  = [Dm2(l  + 32)2  + (GJf  + EIfB2)  + 

(4.4) 

12DZ2  + FV-,  tt2 

ir4m2  A L2  R(l/2  + B2) 

where  pcg(n,m)  is  the  buckling  pressure  for  a specified  n and  m and 
Af  = 1 + 2n2d  (1  - g2y)  + n4d2  (1  + g2)2 

A = (1  + B2)2  + 2B2F  (1  + u)  + B4  F(1  - y2) 

(4.5) 

Z2  = L$  (1  - y2)  / R2t2 
8 = nLs  / mirR 
F = (bw  + bh)  / tL' 

D = Et3  / 12  (1  - y2)  (4.6) 

d = d / R 

where  d is  the  (algebraic)  distance  to  the  neutral  axis  of  the  frame 
from  the  midplane  of  the  hull  plating,  taken  as  positive  when  it  is 
outward  from  the  central  axis  of  the  hull  [25].  The  quantities  GJf 
and  Elf  are  the  torsional  and  bending  stiffness  of  the  frame,  respec- 
tively. Symbols  representing  hull  and  frame  dimensions  are  defined 
in  Figure  9.  The  quantities  G,  E,  and  y are  the  shear  modulus 
G = E / 2(1  + y),  tensile  modulus,  and  poisson's  ratio,  respectively. 


Buckling  of  the  shell  between  frames  is  controlled  by  the  con- 
straint 


cs 


» 0 


where 


g2  - (P*cs  - S2p)  / P* 


(4.7) 
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Pcs(n,m)  = (tt2/L,2R)  [Dm2(l  + 02)2  + l/A(12DZ2/Tr2m2)]l/(l/2  + 02)  (4.8) 

For  a given  set  of  design  parameters  and  pressures  the  buckling 
pressure  functions  p__(n,m)  and  prc(n,m)  depend  on  the  two  discrete 
independent  variables,  the  axial  and  circumferential  wave  number  m 
and  n respectively.  To  obtain  the  critical  buckling  pressure  one  must 
find  the  minimum  of  PCg(n,m)  and  Pcs(n,m)  with  respect  to  n = 0»  1»2.» 

and  m = l,  2 For  this  purpose  the  procedure  of  Reference  [26] 

is  used  here. 

The  necessity  to  avoid  plating  yield  produces  the  constraint 

g3  ’ (apa  - V / °pa  » 0 <4'9> 

where  a is  the  allowable  plating  stress  and 
pa 

CTp  = (cV  + al  ' 0rc<j))1/2  (4J°) 

The  quantities  ar  and  a ^ are  found  from 

= -PR/tn  + T(Hm.+  uHe)]  (4.11) 

°r  = *PR/2t[1  + 2r  HE]  (4.12) 

The  terms  -pR/t  and  -pR/2t  are  the  hoop  and  axial  stresses,  respec- 
tively; r is  a frame  deflection  parameter,  and  and  are  the 
functions  that  define  the  bending  effect  on  the  shell  due  to  local 
frame  reinforcing  [24].  Equations  for  H^,  and  are  given  in  Refe- 
rence [24]. 

The  constraint  against  the  frame  yielding  failure  mode  is  given 


by 


g4  = (afa  " al]  ' afa 


(4.13) 


48. 


* 


| 


i 

I 


where  cr^a  is  the  allowable  frame  stress  and  the  maximum  frame 
stress.  The  quantity  ctt  is  given  by 

°T  = ab  + ac  (4.14) 

with  the  compressive  hoop  stress  ac  given  by 

ac  = Qp  PR/(A  + bt)  (4.15) 

Here  A is  the  cross-sectional  area  of  the  frame  and  Qp  is  the  total 
radial  load  per  inch  of  circumference.  Qpis  given  by  equation  (23) 
of  Reference  [24].  The  frame  bending  stress  ab  results  from  slight 
out-of-roundness  of  the  hull.  It  is  calculated  using 

= [Ece  (n2  - 1)  / R2]  [p/  (pc$  - p)]  (4.16) 

where  c is  the  distance  from  the  midplane  of  the  shell  to  the  surface 
of  the  frame,  e is  the  eccentricity  from  the  true  circle  radius,  and 
n the  wave  number  that  minimizes  p„„. 

It  should  be  noted  that  ab  is  discontinuous  at  points  near  p = pcg 
and  is  discontinuous  with  respect  to  n,  since  n is  an  integer. 

Flange  buckling  is  controlled  by  letting 
95  = ob  - [.5tt2E/12  (1  - v2)]  [x3/  (x4  - x2)]2  * o (4.17) 

Two  geometric  constraint  equations  can  be  used.  A maximum  flange 
thickness  is  specified  to  prevent  the  flange  from  becoming  excessively 
thick.  It  is  desirable  to  specify  a maximum  flange  thickness  as  a 
fraction  of  the  plating  thickness.  Thus,  the  constraint 

g6  = (C2t  - b1)  / C2t  >.  o (4.18) 

is  used,  where  C2  is  the  maximum  flange  to  plating  thickness  ratio. 
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The  constant  C?  is  made  arbitrarily  large  if  the  designer  does  not 
wish  to  apply  this  constraint.  Minimum  flange  width  is  controlled  by 
the  inequality 

g7  = (w  - C3h)  / w ^ o (4.19) 

to  insure  that  the  flange  width  is  sufficient  so  that  the  design  rule 
for  the  control  of  the  web  cribbing  is  not  invalidated,  C3  is  an 
arbitrary  constant  controlling  the  minimum  ratio  of  flange  width  to 
web  height. 


It  is  desirable  to  limit  the  minimum  or  maximum  (depending  on 

whether  internal  or  external  frames  are  used)  radius  of  the  frame 

because  of  space  consideration.  One  can  then  use  the  constraint 

g8  = (R  ♦ h ♦ b'  ♦ t/2  - R^)  / Rm1„  » 0 (4.20) 

where  R.  is  the  minimum  specified  radius  for  internal  frames  or 
mm 

(Rnax-R-t/2-h-b'>  ' W°  <4'21> 

for  external  frames  where  is  the  maximum  specified  radius. 

max 

To  insure  against  web  buckling  the  following  design  constraint 
is  invoked. 

g9  = ab  - [4it2E/12  (1  - v2)]  (x2  / x6)2  ^ 0 (4.22) 


Side  constraints  of  the  form 


(x“  - xi)  / x!|  * 0 


(xi  - x![)  / xi  ^ Oj 

are  used  to  limit  the  range  of  the  variables  x;.  for  manufacturing, 
space,  or  other  practical  reasons.  The  flange  width  clearly  cannot 
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4.3.1  Code  description.  A FORTRAN  IV  code  called  CAD0P4  based 
on  this  algorithm  was  developed  and  used  in  this  study.  The  code  is 
operational  with  IBM  FORTRAN  levels  G and  H and  the  UNIVAC  TDOS  sys- 
tems. The  user  is  required  to  supply  the  objective  and  constraint 
functions,  the  initial  values  of  the  variables,  and  side  constraint 
values. 

The  initial  step  size  is  internally  generated  so  that  at  the 

starting  point  a step  of  size  a0  in  the  Vf  direction  produces  a one 

percent  change  in  the  value  of  the  objective  function.  This  step  is 

constrained  so  that  a0  > .01.  The  minimum  step  size  is  defined  as 

“min  = a°/T000.  The  program  uses  e^  = 10"5,  e2  = e3  = 10”®,  = 6, 

5 

N2  = 10,  = 10  and  the  activity  limit  constant  is  initially  speci- 

fied as  Kj  * 0.01  for  all  constraints. 


exceed  the  frame  spacing,  imposing  an  upper  limit  on  x^.  Sometimes 
x^  is  limited  by  design  consideration,  thus  if  x^j  is  the  designer 
specified  maximum 


The  formulation  of  the  example  treated  here  is  now  complete,  but 
other  forms  of  constraints  and  objective  functions,  which  do  not 
arise  in  this  example  are  clearly  possible. 

4.3  Discussion  of  Results 
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Nondimensional  constraint  equations  of  the  form 

9j (xi ) = (Uk  - Bk)  / Uk  > 0 when  f 0 and  Bk  f 0 
or  (4.25) 

9j(xi ) = (Bk  - Lk)  / Lk  > 0 when  Lk  f o and  Bk  f-  0 
are  used.  Otherwise,  a dimensional  form  of  equations  is  used. 

A constraint  given  as  0 $ b(x.j)  $ A would  be  written  as  two  con- 

\ 

straint  equations.  For  example,  one-could  be  of  the  form  of  the  first 
of  equations  (4.25)  where  B-j  = b(x^)  and  U-j  = A.  The  second  would  be 
of  the  form  gg  = -b(xi)  since  Lk  = 0. 

4.3.2  Code  application.  The  above  problem  was  treated  using 
CAD0P3  and  CAD0P4  codes.  Except  that  the  initial  step  size  is  speci- 
fied for  this  problem  as  a = 4/R,  where  R is  the  radius  of  the  shell. 
The  minimum  step  size  is  defined  as  = a^/200.  Runs  were  made  on 

an  IBM  370  model  168  using  FORTRAN  level  G.  The  same  starting  points 
were  used  in  these  runs.  The  execution  time  in  seconds  along  with 
other  necessary  information  were  printed  at  each  base  point  in  order 
to  closely  track  the  synthesis  path  of  each  procedure. 

Tables  1 and  2 present  the  results  obtained  from  CAD0P3  and 
CAD0P4  codes,  respectively.  A total  of  72810  function  evaluations 
were  performed  using  CAD0P3  code  prior  to  termination.  Due  to  the 
great  execution  time  required  by  CAD0P3  for  this  problem,  the  run  was 
terminated  prior  to  reaching  the  optimum.  The  execution  time  for 
this  run  was  56.2  seconds.  CAD0P4  code  required  only  about  4000  func- 
tion e/aluations  with  4.8  seconds  execution  time  to  reach  a similar 


level  of  convergence. 


The  constraint  presenting  the  web  buckling  (g^)  was  active  from 
near  the  beginning  of  the  search,  while  the  plate  yielding  (g^)  be- 
came active  after  a relatively  short  time  and  remained  active  through 
the  rest  of  the  search.  The  general  buckling  constraint  was  also 
active  at  the  latter  part  of  the  search.  Thus,  the  optimal  design  is 
constrained  by  general  instability,  local  web  buckling  and  plating 
yield. 

The  CAD0P3  code  came  within  1%  of  the  optimum  point  after  23  sec- 
onds and  approximately  43000  function  evaluations,  while  CAD0P4  re- 
quired only  1.4  seconds  and  about  1000  function  evaluations  to  reach 
this  level  of  convergence.  It  should  be  pointed  out  that  for  most 
practical  engineering  purposes,  a 1%  level  of  accuracy  is  quite  accep- 
table. It  is  the  refining  process  that  is  responsible  for  the  greatest 
portion  of  the  total  execution  time. 

It  must  be  pointed  out  that  due  to  the  presence  of  intermediate 
print  statements,  the  execution  time  required  for  solution  without 
these  statements  is  appreciably^ ess  than  the  above  mentioned  time. 

A run  made  without  any  intermediate  print  statements  using  CAD0P4 
code  terminated  after  2.1  seconds  at  the  optimum  point  while  the 
original  run  required  8.7  seconds. 

Comparison  of  the  above  results  demonstrate  the  apparent  supe- 
riority of  the  Boundary  Tracking  algorithm  to  the  DSFD  procedure  with 
respect  to  speed.  This  is  the  result  of  great  reduction  in  the  number 
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of  function  evaluations  required  for  convergence.  It  demonstrates 
the  ability  of  the  BT  procedure  to  move  efficiently  along  a difficult 
behavior  boundary. 


J 


CHAPTER  5 


COMPARISON  STUDY 

The  CAD0P4  code  as  described  in  section  4.3.1  was  used  in  this 
study  along  with  all  the  constants  and  parameters  presented  earlier. 

This  comparison  study  is  based  on  the  works  of  Eason  and  Fenton 
given  in  References  [14]  and  [15].  All  ten  problems  treated  in  these 
references  were  run,  using  the  new  code  with  the  starting  points 
given  in  Reference  [14].  All  the  control  constants  and  the  step 
sizes  are  internally  generated  so  that  there  was  no  special  tuning 
for  individual  problems.  Computations  were  performed  in  double  pre- 
cision on  an  IBM  370  model  168  system,  using  a level  G compiler,  thus 
closely  simulating  the  Eason  and  Fenton  study. 

A brief  description  of  the  codes  studied  are  given  in  Table  3, 
while  Table  4 contains  the  data  for  rating  the  success  of  the  various 
codes  in  solving  the  ten  problems  to  which  they  were  applied.  A de- 
tailed description  of  the  problems  is  given  in  Reference  [14].  A 
normalized  time  required  for  the  solution  of  each  problem  successfully 
solved  is  given  in  Table  4.  Normalized  time  is  defined  here  as  the 
execution  CPU  time  divided  by  the  CPU  time  required  to  execute  a tim- 
ing standardization  program  [14].  The  symbol  "P"  denotes  progress 
toward  a solution,  and  a blank  denotes  failure.  The  criteria  used  in 
References  [14]  and  [15]  for  defining  successful  solution  or  progress 
were  also  applied  here. 


r 

1 

t 

• 

TABLE  3 
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CODE  DESCRIPTIONS 

Name 

Description 

Algor 

Class* 

ADRANS 

Random  search  followed  by  pattern  moves 

DS 

CLIMB 

Rosenbrock  search 

DS 

DAVID 

Davidon-Fletcher-Powell  with  numerical  derivatives 

GF 

DFMCG 

Fletcher-Reeves  conjugate  gradient  method  with 
secant  approximation  derivatives 

G 

DFMFP 

Davidon-Fletcher-Powell  with  secant  approsimation 
derivatives 

ii 

G 

FMIND 

Hook  & Jeeves  pattern  search 

DS 

GRAD4 

Steepest  descent  method 

G 

GRID1 

Grid  and  star  network  search,  with  shrinkage 

AR 

MtMGRD 

Davidon-Fletcher-Powell  with  retained  step  length 
information 

GF 

NMSERS 

Simplex  search 

DS 

PATSH 

Modified  pattern  search,  dome  strategy 

• DS 

PATRNI 

Modified  pattern  search,  ridge  strategy 

DS 

RANDOM 

Random  search  with  shrinkage 

AR 

SEEK1 

Pattern  search  followed  by  random  search 

DS 

SEEK3 

Modified  pattern  search 

DSF 

SIMPLX 

Modified  simplex  search 

DSF 

DSDA 

Modified  pattern  search  followed  by  Mugele's  search 

DS 

CAD0P2 

Modified  pattern  search  followed  by  Zoutendijk 
feasible  direction  method 

DS  | 

CAD0P3 

Modified  CAD0P2 

DS 

CADOP4 

1 , 

Boundary  Tracking  method 

BT 

*DS  = direct  search,  DSF  = direct  search  e ploying  SUMT  strategy  and 
penalty  function,  G = gradient  procedure,  GF  = gradient  procedure 
using  SUMT  strategy  and  penalty  function,  AR  = area  reduction  method, 
BT  3 Boundary  Tracking. 


TABLE  4 

PERFORMANCE  OF  OPTIMIZATION  CODES* 


regional  constraints. 


The  Important  variables  effecting  execution  time  are;  the  prob- 
lem, the  algorithm,  the  code,  the  input  and  output  requirements,  the 
architecture  of  the  machine  and  the  compiler  system.  In  order  to 
minimize  the  effects  of  those  variables  not  involved  in  the  compari- 
son, similar  computer  machines  and  compiler  systems  along  with  the 
same  test  problems  and  output  requirements  are  used  here.  Unfortuna- 
tely the  effects  of  the  code  structure  can  not  accurately  be  deter- 
mined. Thus,  data  of  Table  4 and  ratings  shown  in  Table  5 can  not 
be  considered  to  accurately  reflect  the  effectiveness  of  the  basic 
optimization  algorithms. 

The  data  and  rating  procedures  used  here  however,  can  be  direct- 
ly applied,  with  reasonable  accuracy,  for  relative  comparison  of  the 
codes  tested.  The  results  obtained  are  thus  useful  to  the  user  of 
such  codes.  The  relative  effectiveness  of  basic  optimization  algorithms 
can  only  be  infered  from  this  data  if  one  assumes  that  the  codes  for 
each  procedure  have  been  prepared  with  essentially  comparable  effec- 
tiveness. 

This  study  uses  Eason  and  Fenton's  comparison  procedures  since 
they  are  the  best  available  and  the  most  current.  Furthermore,  use 
of  another  investigator's  comparison  methods  rather  than  one  proposed 
by  the  developer  of  a new  method  reduces  the  tendency  toward  developer 
bias  in  reporting  on  the  effectiveness  of  his  method. 

The  data  from  Table  4 may  be  applied  to  a number  of  rating  methods 
for  comparing  the  codes  tested.  Table  5 presents  relative  ranking  of 


RELATIVE  RANKING  OF  OPTIMIZATION  CODES 


CLIMB  63.6  DAVID 

DFMCG  92.5  ADRANS 
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the  codes  using  the  criteria  of  Reference  [15].  The  rating  equations 
are  as  follows: 


Let  the  number  of  problems  solved  by  code  "a"  be  n , and  let  n1 

a d 

be  the  number  of  problems  with  a "P"  rating.  Then  the  numerical  suc- 


cess rating  N,  is  given  by 

a 


Na  = "a  + 


(5.1) 


One  of  two  computational  efficiency  ratings  is  given  as 

fa  ' [£J=1  bap  bap  ' <‘ap>]  ' "a  <5-2> 

where  b = 1 if  code  "a"  solved  problem  "p"  and  zero  otherwise,  t 

dp 

is  the  normalized  time  required  for  solution  and  mi n ( tap ) is  the 
shortest  time  required  by  any  of  the  codes  studied  to  solve  problem 


'p".  The  other  efficiency  criterion  is 


fa  = [£p=l  bap  bap  1 raean  (tap):l  1 \ 


(5.3) 


where  mean(tap  is  the  average  time  required  by  the  codes  studied  to 
solve  problem  "p".  An  overall  rating  number  which  can  be  considered  a 
composite  measure  of  generality  (reliability)  and  efficiency  (speed) 


is  given  by 


T = fi 
‘a  Vl  ap 


(5.4) 


where  tap  is  set  equal  to  t if  algorithm  "a"  solves  a problem  "p", 
and  to  twice  the  time  used  by  the  slowest  code  solving  a problem  "p" 
if  code  "a"  could  not  solve  it.  This  penalty  time  is  used  to  penalize 
code  unreliability.  Only  codes  that  solved  half  or  more  of  the  prob- 
lems are  rated  for  efficiency. 

The  tables  presented  here  are  similar  to  those  given  in  References 
[15]  and  [17],  except  that  the  CAD0P3  and  CAD0P4  codes  are  included. 
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The  CAD0P3  code  is  a modified  version  of  CAD0P2  code  presented  in  Ref- 
erence [17].  Thus,  the  rating  methods  and  presentation  of  results  used 
are  essentially  those  of  Reference  [15]  and  [17]. 

From  Tables  4 and  5 it  can  be  concluded  that  CAD0P4  is  the  fastest 
code  in  overall  generality  and  efficiency  rating.  This  Boundary  Tracking 
method  presented  here  appears  comparable  in  speed  to  the  fastest  methods 
presented  in  Reference  [17],  (NMSERS,  PATRNI,  DSDA,  and  CAD0P2).  CAD0P4 
is  faster  than  average  in  all  problems.  It  is  the  fastest  code  solving 
problems  1,  2,  3,  6,  7,  and  9,  which  are  problems  with  behavior  con- 
straints. On  the  other  hand,  problems  4,  5,  8,  and  10,  where  the  CAD0P3 
and  CAD0P4  codes  performed  identically,  but  less  efficiently  than  some 
other  procedures,  are  problems  without  such  constraints.  This  is  ex- 
pected since  the  two  codes  are  identical™  treating  the  unconstrained 
problems. 

Of  the  problems  tested,  problem  1 which  has  the  largest  number  of 
behavior  constraints  (10)  was  apparently  also  the  most  difficult,  since 
it  was  solved  by  only  six  of  the  twenty-one  codes  presented  in  Table  4. 
The  second  most  difficult  problem  was  problem  9,  for  which  only  7 of  the 
21  codes  led  to  the  optimum  point.  This  problem  has  only  6 behavior 
constraints  and  two  variables.  However,  there  is  a rapid  change  in  the 
slope  of  the  behavior  boundary  in  the  region  near  the  optimum  point. 
Because  of  this  rapid  variation  in  slope,  a major  part  of  the  total  ex- 
ecution time  was  spent  on  refining  the  location  of  the  point  so  as  to 
meet  the  convergence  criterion.  Problem  10,  which  required  the  longest 
execution  time,  has  no  behavior  constraints.  The  long  execution  time 
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required  by  this  problem  ts  due  to  the  complexity  of  the  objective 
function  and  the  number  of  the  function  evaluations  needed  to  reach 
the  optimum. 

It  can  be  seen  from  Table  4 that  the  CPU  time  required  for  trea- 
ting problems  9 and  10  in  most  cases  is  much  greater  than  the  sum  of 
the  CPU  times  required  for  the  solution  of  the  rest  of  the  problems 
tested.  Thus,  the  speed  of  a code  in  solving  problems  9 and  10  plays 
an  overwhelming  role  in  determining  the  code's  overall  generality  and 
efficiency  rating,  T . For  example,  a code  that  is  very  efficient  on 
the  majority  of  problems  but  requires  long  solution  times  for  problems 
9 and  10  would  have  a relatively  low  overall  generality  and  efficiency 
rating  using  equation  (5.4).  Therefore,  this  rating  method  is  not 
very  useful  for  comparing  codes  solving  all  the  test  problems.  For 
such  codes  the  efficiency  ratings  f,  and  f,  represent  the  ability  of 

a a 

the  method  more  clearly.  Furthermore,  the  overall  generality  and 
efficiency  rating  T,  does  not  sufficiently  penalize  those  codes  which 

a 

failed  to  solve  any  of  the  test  problems  1 through  8,  but  performed 
well  on  problems  9 and  10.  The  penalty  time  used  here  fails  to  accu- 
rately take  into  account  the  weaknesses  of  such  codes  in  solving 
problems  1 through  8. 

The  superior  performance  of  CAD0P4  code  in  problems  with  behavior 
constraints  compared  to  CAD0P3  code  is  due  to  the  fact  that  the  number 
of  function  evaluations  required  in  CAD0P3  is  proportional  to  the  num- 
ber of  bases  needed  for  solution  multiplied  by  the  number  of  variables. 
On  the  other  hand,  the  number  of  function  evaluations  in  CAD0P4  is 
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proportional  only  to  the  number  of  base  points,  and  does  not  have  to 
be  multiplied  by  the  number  of  variables.  Therefore,  in  CAD0P3  the 
number  of  function  evaluations  increases  directly  with  the  number  of 
variables,  while  in  the  CAD0P4  code  it  does  not. 

Compared  to  the  relatively  reliable  codes  (n.  > 8),  CAD0P4  is 

a 

faster  than  CAD0P3  on  six  of  the  ten  problems  solved  by  CAD0P3,  faster 
than  DSDA  on  six  of  nine,  and  faster  than  PATSH  on  eight  of  nine,  and 
faster  than  SEEK3  in  all  problems  solved  by  these  schemes.  As  may  be 
seen,  the  CAD0P2  and  CAD0P3  codes  are  the  only  reliable  codes  compar- 
able to  CAD0P4.  CAD0P4  is  significantly  faster  than  CAD0P3,  CAD0P2, 
PATSH,  SEEK3,  and  SIMPLX.  The  straight  pattern  (PATRN1 ) or  simple 
(NMSERS)  codes  are  ranked  relatively  high  in  efficiency  primarily  due 
to  their  superior  performance  in  unconstrained  problems  (Table  5). 

Of  this  group,  only  NMSERS  appears  to  be  sufficiently  reliable  to  merit 
consideration  for  use.  The  minor  difference  in  efficiency  between 
NMSERS  and  CAD0P4  is,  however,  overshadowed  by  the  superior  reliability 
of  the  latter.  Thus,  in  the  overall  speed,  generality,  and  efficiency 
rating,  CA00P4  stands  out.  Viewed  on  the  basis  of  these  comparisons, 
CAD0P4  appears  to  be  a superior  nonlinear  mathematical  programing  code. 

The  new  algorithm  is  intended  to  be  a constrained  optimization 
algorithm.  Therefore,  the  performance  of  CAD0P4  is  compared  to  the 
codes  of  References  [15]  and  [17]  on  problems  having  behavior  con- 
straints. This  comparison  is  shown  in  Table  6. 

Table  6 is  similar  to  Table  5 except  that  the  problems  without 
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behavior  constraints  are  excluded.  The  rating  equations  used  for  these 
values  are  the  same  as  those  used  for  Table  5.  The  results  of  Table  6 
indicate  the  superiority  of  CAD0P4  with  respect  to  speed,  generality, 
and  efficiency.  Only  codes  that  solved  half  or  more  (n.  ^ 3)  of  the 

d 

problems  are  rated  for  efficiency. 

The  superior  performance  of  the  CAD0P4  code  in  all  the  above  ra- 
tings strongly  suggests  the  superiority  of  the  BT  method  for  treating 
problems  with  behavior  constraints.  CAD0P4  proved  to  be  substantially 
more  efficient  than  the  CAD0P3  and  NMSERS  codes  on  such  problems.  In 
the  overall  generality  and  efficiency  rating  CAD0P4  code  again  stands 
alone  with  CAD0P3  code  in  the  second  place  being  approximately  three 
times  slower.  However,  as  explained  earlier,  due  to  the  relatively 
long  execution  time  required  for  problem  9 this  rating  method  does  not 
reflect  accurately  the  relative  effectiveness  of  the  codes  which 
solved  all  the  test  problems.  The  efficiency  ratings  f.  and  T may 

a a 

again  be  more  useful  in  comparing  such  codes.  In  these  ratings  CAD0P4 
code  is  two  to  three  times  faster  than  other  relatively  fast  codes 
(CAD0P3,  NMSERS). 

Thus,  the  Boundary  Tracking  algorithm  presented  here  appears  to 
be  superior  to  DSFD  with  regard  to  efficiency  and  to  all  other  optimi- 
zation procedures  tested  with  respect  to  generality  and  efficiency. 

Based  on  its  performance  on  the  ten  problems  of  Reference  [15] 
and  on  the  other  comparison  studies,  the  Boundary  Tracking  method 
appears  to  be  a fast  and  reliable  nonlinear  mathematical  programing 
optimization  procedure. 
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In  conclusion,  it  must  be  pointed  out  that  the  test  problems  used 
in  the  above  comparison  study  are  relatively  small  and  simple,  while 
many  practical  engineering  problems  are  large  and  contain  complex  and 
computationally  demanding  functions.  Therefore,  although  it  is  reason- 
able to  assume  that  the  above  comparison  study  is  representative  of 
the  performance  of  the  new  code  on  most  practical  problems  of  size 
and  complexity  similar  to  that  of  the  study,  this  performance  may  not 
be  representative  for  large  complex  problems.  Unfortunately,  no 
comparison  study  utilizing  large  complex  problems  is  available,  and 
is  not  likely  to  be  available  soon  due  to  the  difficulty  and  cost 
associated  with  such  a project.  Furthermore,  one  cannot  guarantee 
that  the  performances  noted  in  the  present  study  are  typical  for  all 
relatively  small  simple  problems.  Nevertheless,  this  is  the  most  com- 
plete of  the  available  general  comparison  studies  and  should  be  useful 
as  a guide  in  selecting  a suitable  algorithm  for  a particular  problem. 
The  designer  must,  however,  proceed  with  care  in  comparing  his  problem 
to  the  above  test  problems,  taking  into  account  such  things  as  the 
number  of  variables,  constraints,  nature  of  the  function  to  be  evalua- 
ted, etc.  in  selecting  a desirable  algorithm. 
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CONCLUSION 

The  successful  application  of  the  CAD0P4  code  to  the  relatively 
complex  problem  of  shell  synthesis  and  its  performance  in  the  general 
comparison  study  presented  in  section  4.1,  imply  that  the  Boundary 
Tracking  method  developed  here  is  a superior  nonlinear  mathematical 
programming  procedure.  Although  the  code  proved  to  be  the  best  in  the 
general  comparison  study,  the  real  potentials  of  the  algorithm  are 
demonstrated  in  the  shell  design  problem.  Since  many  engineering 
problems  are  of  such  a class,  the  new  algorithm  shows  the  promise  of 
adding  a major  contribution  to  the  field  of  automated  design. 

In  the  above  studies  tl)e  BT  method  proved  to  have  great  potential 
for  use  in  computationally  demanding  problems.  However,  as  in  the  case 
of  most  new  methods,  additional  studies  may  lead  to  improvement,  par- 
ticularly in  the  speed  of  the  algorithm. 

The  only  apparent  disadvantage  of  CAD0P4  is  its  relative  complex- 
ity, compared  to  some  of  the  other  reasonably  reliable  methods.  It 
contains  1320  FORTRAN  statements,  while,  for  example,  the  CAD0P3  con- 
tains 790,  DSDA  contains  372,  and  PATSH  only  75  FORTRAN  statements. 
Thus,  if  CAD0P4  code  is  not  available  in  a compiled  form,  for  simple 
problems  of  relatively  low  complexity  and  dimensionality,  one  of  the 
simpler  codes  may  be  preferable  since  the  time  required  for  compilation 
of  the  program  may  exceed  the  time  saved  by  using  CAD0P4. 
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In  addition,  its  performance  on  problems  without  behavior  con- 
straints, although  quite  good,  is  not  as  outstanding  as  on  problems 
with  such  constraints. 

In  conclusion,  the  user  of  this  or  any  other  algorithm  for  auto- 
mated design,  should  be  aware  that  the  general  nonlinear  constrained 
optimization  problem  is  quite  difficult  to  handle.  Also,  none  of  the 
available  techniques  will  guarantee  an  optimum  solution.  One  should, 
therefore,  be  careful  in  the  application  of  such  an  algorithm  and 
analyze  the  results  thoroughly.  One  also,  should  make  use  of  several 
synthesis  runs,  using  different  starting  points  where  possible,  before 
assuming  the  value  is  an  optimum  solution. 
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USER  INSTRUCTIONS  FOR  CADOP4 
I.  INSTRUCTION 

This  code  treats  inequality  constrained  or  unconstrained 
linear  or  nonlinear  minimization  problems  by  means  of  a direct 
search  mathematical  programming  method.  The  method  starts  from 
a user  specified  starting  or  initial  point  and  generates  a se- 
quence of  better  points  until,  hopefully,  an  optimal,  or  near 
optimal,  point  is  reached.  The  code  is  intended  primarily  for 
constrained  nonlinear  problems.  Linear  problems  are  much  more 
effectively  handled  by  one  of  many  linear  programming  methods. 
Unconstrained  problems  are  better  treated  by  UNCDP,  a simplified 
version  of  CAD0P2,  or  one  of  several  other  unconstrained  non- 
linear problem  codes.  CAD0P4  will,  however,  treat  both  linear 
and  unconstrained  problems. 

It  should  be  noted  that  nonlinear  mathematical  programming 
methods  generally  do  not  guarantee  an  optimal  solution  because 
of  the  possibility  of  multiple  local  optima,  algorithm  failure 
or  numerical  difficulties.  Although  CAD0P4  is  among  the  most 
reliable  of  the  nonlinear  mathematical  programming  codes  it  can- 
not be  considered  absolutely  reliable.  Thus,  it  is  desirable 
to  use  several  optimization  runs  with  widely  separated  starting 
points  to  confirm  the  convergence  has  been  achieved  or  to  deter- 
mine if  local  optima  are  present. 


This  code  treats  the  problem: 

Find  those  values  x^  of  the  set  of  "variables"  x^, 

i = 1,2 IP  and  the  constant  "parameters"  p^»  k = 1,2 KP , 

that  result  in  the  minimum  value  of  the  objective  function 


f (Pk,xi),  that  is 


f ^pk,xi^  = min  k ^P»xi^ 
while  also  satisfying  the  "behavior"  constraints 


gj  (pkXi)  * 0 i “ 1>2>— Jp 

and  the  "regional"  constraints 

1 


"i  xi  ' "i 

The  program  in  its  present  form  treats  problems  with 


x j 


2 < IP  < 10,  0 < JP  < 10  and  0 < KP  < 100 
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III . SPECIFICATION  OF  OBJECTIVE  AND  CONSTRAINT  FUNCTIONS 

The  subprogram  FUNCTION  OBJ(J)  is  essentially  a dummy  sub- 
program created  to  accept  the  users  FORTRAN  IV  program  state- 
ments defining  the  objective  and  behavior  constraint  functions. 
The  normal  form  of  the  behavior  constraints  is 

Bj  ^pk’xi^  - Uj  ^pk,xi^ 

where  can  be  thought  of  as  the  behavior  and  the  upper 
limit  of  behavior. 

The  FUNCTION  OBJ(J)  Is  modified  so  as  to  define  the  objec- 
tive and  constraint  functions  by  inserting  the  following  state- 
ments immediately  after  the  COMMON  statement 

GO  TO  (1,2, jp),  J (omit  if  no  constraints  are  used) 

OBJ  * expression  defining  f (P^*3^)  using  arrays  P(k),  X(i) 
RETURN 

j B * expression  defining  B ^ 

U - expression  defining  U^ 

GO  TO  101 

A constraint  statement  group  is  used  for  every  constraint.  The 
last  constraint  group  need  not  use  the  final  GO  TO  statement. 
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Example 

For  the  problem: 
minimize 


2xl  x2  " P1X3X4‘ 


such  that 


2x3  - x1  <_  p. 


and 


2 

-6  < x0  - x,  < 0 

— 2 4 — 

one  could  define  three  constraints  and  insert  the  statements 
GO  TO  (1,2,3),  J 

OBJ  = 2.  * X (1)  * X (2)  - P (1)  * X (3)  * X (4)  **2 
RETURN 

1 B = 2.*  X (3)  **  3 - X (1) 

U = P ( 2 ) 

GO  TO  101 

2 B - X ( 2 ) **  2 - X(4) 

U - 0. 

GO  TO  101 

3 U - (X  (2)  **  2 - X (4)  ) 

B - - 6. 
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Alternately,  noting  that  both  upper  and  lower  limits  on 
the  second  constraint  equation  cannot  be  active  simultaneously 
one  could  use  two  constraint  equations  and  insert  the  program 
segment : 

GO  TO  (1, 2) , J 

ORJ  = 2.  * X (1)*  X ( 2 ) - P (1)  * X (3)  * X (4)  **  2 
RETURN 

1 B =*  2.  * X (3)  **  3 - X (1) 

U = P ( 2 ) 

GO  TO  101 

2 TB  = X (2)  **  2 - X (4) 

LF (TB . LT . -3 . ) GO  TO  3 

B = TB 

U - 0 

TO  TO  101 
3 B - -6. 

U - TB 

The  contraints  values  at  the  optimum  are  printed  in  the 
form  of  equation  (2)  where 

f_(Bj  ‘ Uj)/|Ujl  Uj  * 0 and  B * 0 
gj  “ 1 

^“(Bj  - Uj ) Uj  = 0 or  Bj  * 0 

thus  a positive  value  of  gj  indicates  the  constraint  is  satisfied. 
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IV  DATA  CARDS 

The  program  allows  multiple  optimization  runs  of  a par- 
ticular problem  within  a single  computer  run  using  different 
parameters,  for  parametric  studies,  different  starting  points, 
for  optimality  confirmation  or  search  for  local  optima,  and  dif- 
ferent sets  of  regional  constraint  limits. 

The  user  must,  by  use  of  the  data  card  set,  specify  the 
number  of  parameters,  variables  and  behavior  constraints  used, 
specify  for  the  first  run  if  regional  constraints  are  to  be 
used,  and  for  subsequent  runs  whether  the  parameters,  initial 
point  or  regional  iimits  are  to  be  changed.  If  parameters  are 
used  they  must  be  entered.  The  initial  point  must  be  entered. 
The  values  for  the  regional  constraints  must  be  entered  if  such 
limits  are  imposed. 

The  first  data  set  is  entered  on  the  "Problem  Constraints" 
card  which  specifies,  in  order,  the  number  of  parameters  (KP), 
variables  (IP)  and  behavior  constraints  (JP),  used  (see  line 
8 of  program  listing) . The  data  is  entered  in  3110  format  in 
the  first  three  fields  in  order  and  must  be  right  .justified. 

The  2nd  data  set  is  entered  on  the  "Data  Control"  card 
and  is  used  to  specify  whether:  1)  new  parameters  are  to  be 

used  (ICNTR2),  2)  a new  initial  point  is  to  be  used  (ICNTR3) 
and,  3)  regional  limits  or  new  regional  limits  are  to  be  used 
(ICNTR4)  (see  lines  12-23).  The  entries  are  made  in  3110  for- 
mat in  the  order  given  above.  For  the  first  run  no  entries  in 


I 


3 


iJ 


77 


| 

: 


the  first  two  fields  are  needed.  If  regional  limits  are  used 
an  entry,  any  nonzero  digit,  is  made  in  the  third  field  (col. 

21- 30)  . 

If,  and  only  if,  the  number  of  parameters  specified  is 
greater  than  zero  (KP  > 0)  a 3rd  data  set  is  used  to  enter  the 
problem  parameters  p^.  The  entries  are  made  5 to  a card  in 
F15.5  format  in  order  i = 1,2  KP , (see  lines  13  and  14). 

The  4th  data  set  specifies  the  starting  point  and  entries 
are  made,  in  order,  8 to  a card  in  F10.5  format  (see  lines  7 
and  15)  . 

i 

If,  and  only  if,  an  entry  is  made  in  the  third  field  of 
the  Data  Control  card  (ICNTR4  ^ 0)  a 5th  data  set  defining  the 
lower  limits  is  entered,  8 to  a card,  in  F10.5  format  followed 
by  an  upper  limit  set  (starting  with  a new  card)  with  similar 
format.  Both  complete  limit  sets  must  be  entered  (see  lines 

22- 25)  . 

For  every  additional  run  an  additional  Data  Control  card 
is  added  followed  by  parameter  and/or  initial  point  and/or 
regional  limit  sets,  where  and  only  where  the  need  for  a new 
set  is  Indicated  by  an  appropriate  entry  (ICNTR2  + 0 or 
(ICNTR3  + 0 or  ICNTRL4  f 0)  on  the  Data  Control  card  (see  lines 
18-22)  . 

A blank  Data  Control  card  is  used  to  terminate  the  run(s) 
(see  lines  18  and  19)  . 
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V.  CONVERSION  TO  SINGLE  PRECISION 


The  program  as  supplied  is  a double  precision  form.  Ex- 
perience has  shown,  however,  that  the  great  majority  of  prob- 
lems can  be  treated  adequately  in  single  precision  and  thus 
it  is  recommended  that  single  precision  be  tried  first.  To 
convert  to  single  precision,  remove  the  IMPLICIT  statement  at 
the  beginning  of  MAIN  and  all  subprograms.  Insert  function 
conversion  statements  such  as  DSQRT(B)  = SQRT(B),  etc.  after 
the  COMMON  block  in  the  subprograms  where  such  functions  are 
used.  Furthermore,  the  REAL  FUNCTION  0BJ*8(J)  should  be  put 
in  single  precision  form. 


Storage  requirements  for  the  basic  program  in  double  pre- 
cision is  about  150K  and  in  single  precision  about  100K. 
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VI.  CHANGE  OF  PROBLEM  SIZE 


To  change  the  maximum  number  of  variables  (IP)  or  the  max- 
imum number  of  constraints  (JP)  the  program  can  treat,  change 
the  arrays  in  MAIN  and  all  subprograms  as  follows : 

1.  Change  arrays  of  the  COMMON  statements  of  the  main  and 
all  subprograms  as  follows; 

D,  DL,  DU,  and  DZ  become  D(IP),  DL(IP),  DU(IP)  and  DZ(IP) 
XB  becomes  XB(IP,5). 

V (A  In  0PT2 , CONMOV  and  LINPRO)  becomes  VCM,N).  Where 

M = JP  + 4lP  + 5 
N - JP  + 5IP  + 6 

GN  (GO  in  0PT2,  CONMOV,  ZERO,  and  G in  PATTRN  and  FIND) 
becomes  GN(JP).  IA  becomes  IA(JP) . 

where  IP  is  the  new  maximum  number  of  variables  and  JP 
the  new  number  of  constraints 

2.  In  COMMON/DOL/  change  BL  array  to  BL(JP)  wherever  this 
statement  is  used. 

3.  In  COMMON/ZE/  change  DT  and  DTT  to  DT(IP)  and  DTT(IP). 
change  SE  to  SE(JP). 

change  DK  to  DK(JP  + 1,  IP). 

4.  In  COMMON/POL/  change  GB  to  GB(JP). 

5.  In  COMMON/CON/  change  SS  to  SS(IP). 

CHANGE  IJ  and  JI  to  IJ(JP)  and  JI(JP). 

6.  In  COMMON/FIN/  change  x to  x(IP). 


change  IA,  CK  and  IC  to  IA(JP) , CK(JP)  and  IC(JP) 
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7.  In  the  SUBROUTINE  DELF  Dimension  statement  change  DT  and 
DTT  to  DT(IP),  DTT(IP)  and  GO  to  GO(JP). 

8.  In  the  SUBROUTINE  DB  DIMENSION  statement  change  DT  to  DT(IP). 

9.  In  the  0PT2  SUBROUTINE  DIMENSION  statement  change  IB,  IC,  and 
X to  IB (IP)  etc..  TN  and  GG  to  TN(JP),  GG(JP).  change  SSS  to 
SSS(M),  where  M is  defined  in  item  1 as  DB1  to  DB1  (JP  + 1, 
IP)  . 

10.  In  the  SUBROUTINE  CONMOV  DIMENSION  statement  change  GG  to 
GG  (JP)  and  X to  X (IP) 

11.  In  the  SUBROUTINE  ZERO  DIMENSION  statement  change  GG  to 
GG  (JP). 

12.  In  the  SUBROUTINE  FIND  DIMENSION  statement  change  DT(10)  to 
DT (IP)  . 

13.  In  the  SUBROUTINE  PATTERN  DIMENSION  statement  change  all 
10's  to  IP. 


14.  Change  all  dummy  arrays  such  as  DUM,  IDUM,  IDUN,  etc.  in  all 
COMMON  statements  to  equalize  COMMON  size. 


C001  1 MP  LICIT  RfcAL»»  (A-H.b-Z ) 
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CVtbS.bb  l,Q,IP,KP, JP.1DUN122) 
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